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Abstract. We consider a nonequilibrium process on a timeline with discrete sites 
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flights in the limit of zero space dimensions. The critical properties are investigated 
by extensive numerical simulations and compared with field-theoretic predictions. 
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1. Introduction 



In the present paper we consider a class of probabilistic spreading processes on a timeline. 
As sketched in Fig. [1], the timeline consists of discrete sites t G Z which can be either 
active or inactive. After specifying a certain initial configuration of active and inactive 
sites along this timehne, the process evolves dynamically by subsequent updates of the 
lattice site according to the following probabilistic rules: 

• If the updated site at time tu is active, it attempts to active n{tu) lattice sites in 
the future, where n{tu) =0,1,2,... is randomly selected from a given distribution 
with a finite average n. 

• For each of these attempts a random time interval At = 1,2, .. . is drawn from 
another probability distribution P{At). If the target site at time tu + At is still 
inactive it will be activated, otherwise nothing happens. 

In the following we are interested in the special case where the time intervals are 
asymptotically distributed by a power law 

p(At) ~ (A^)-l"^ (1) 

where < k < 1 is a control parameter. In this case the process can be considered as a 
model for the spreading of activity on a timeline by means of temporal Levy flights [1,2]. 

Remarkably, the timeline spreading process (TSP) shown in Fig. [T] displays a 
dynamical non-equilibrium phase transition when the parameter n is varied. For 
example, if n is sufficiently small the process starting with a single active site at t = will 
terminate after some time, while for large n the process may survive forever, producing 
an asymptotically constant density of active sites. As will be shown below, the two 
regimes of survival end extinction are separated by a well-defined critical threshold Uc, 
where the system undergoes a continuous phase transition. Despite the simplicity of the 
model this transition turns out to be characterized by a surprisingly non-trivial critical 
behavior [3,4]. 

The other interesting aspect of the TSP is the non-Markovian nature of its 
dynamics. Since activity may spread over arbitrarily long time intervals At, the actual 
state of an updated site depends not only on the previous time step but rather on the 
entire history of the process. Therefore, the initial condition is not determined by the 
state of the system at t = alone, instead one has to specify the complete configuration 
of active sites along the entire timeline. In fact, the TSP is probably the simplest model 
which allows one to study non-Markovian features in the context of continuous phase 
transitions far from equilibrium. 

The transition in the TSP belongs to the category of so-called absorbing phase 
transitions [5-9] because it may be interpreted as a transition from a fiuctuating 
active state into a frozen inactive state where the process terminates. Such transitions 
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Figure 1. Timeline spreading process (TSP): Example of a temporal evolution starting 
with a single active site at < = 0. Active sites are updated from left to right, as marked 
by the green arrows. The blue arrows indicate how activity spreads by Levy flights 
over randomly chosen displacements At. 



are associated with certain universality classes, the most prominent one being the 
universality class of directed percolation (DP) [10], which plays a similar role as the 
Ising class in equilibrium statistical mechanics. It is represented e.g. by the contact 
process [11] which is a toy model for epidemics where activity spreads to nearest 
neighbors on a dimensional lattice by means of a Markovian update rule. Although 
this model is easy to define it could be realized experimentally only one year ago by 
Takeuchi et al [12]. Recently DP was generalized to include long-range interactions by 
incorporating spatial [13-15] and temporal [16] Levy flights as well as a combination of 
both [17] (for a review see [18]). As we will see, the TSP studied here may be interpreted 
as the zero- dimensional limit of DP with temporal Levy flights because here we have a 
single site that evolves in time. 

From a broader perspective, the TSP is also of conceptual interest. According to 
a well-known theorem by Landau, phase transitions in equilibrium models with short- 
range interactions require at least two space dimensions. In the non-equilibrium case, 
however, phase transitions in models with short-range couplings such as DP are possible 
in one spatial dimension. The TSP demonstrates that by introducing a non-Markovian 
(i.e. temporally non-local) dynamics phase transitions are possible even in zero space 
dimensions. 
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Figure 2. Typical temporal evolution of the TSP for k = 0.5 on a logarithmic time 
scale. The red dots indicate active sites while the half circles illustrate the Levy flights 
directed forward in time. As can be seen, activity occurs in form of intermittent bursts. 

The paper is organized as follows. In the following section we introduce the specific 
variant of the TSP studied in this paper and discuss phenomenological properties of the 
phase transition. In Sect. 3-4 we summarize previous field-theoretic results [3] which al- 
low one to identify a mean-field and a fluctuation-dominated regime and to compute the 
critical exponents to one-loop order. In Sect. 5 we confirm these findings by numerical 
simulations whereas the conjugate field and the special role of the survival probabil- 
ity will be discussed in Sect. 6. The relation to a recently studied boundary-induced 
phase transition in 1+1 dimensions [3,4] will be discussed in a forthcomming publication. 



2. The model and its phenomenological properties 
2.1. Definition and numerical implementation 

There are many possible variants of the TSP which differ mainly in their specific 
probability distribution for the number of attempted activations n{t) and in the short- 
range details of the Levy distribution P[At). These variants may have different critical 
thresholds but for given k G (0, 1) their critical properties at the transition are 
expected to be universal. 

The numerical results reported in this paper were obtained by simulating the 
following variant of the TSP. It is defined on a timeline with discrete sites s{t) at time 
t G Z, where the values s{t) = 0, 1 denote inactive and active sites, respectively. For 
each run one first has to specify the initial configuration by assigning certain values to 
all sites. In most cases we will start with a single seed of activity at t = by setting 
s{t) := 6tfi. After specifying the initial state an update loop over all time steps t^ is 
executed which starts with the minimal t for which s(t) = 1 and ends when the process 
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terminates, limited by some cutoff time tmax- Inside this loop the following steps are 
carried out: 

(i) If s{tu) = go to (v). 

(ii) Generate a random number 2; G (0, 1). If 2 > p go to (v). 

(iii) Generate another random number y E (0, 1) and set At := y~^/'^. 

(iv) If t„ + At < tmax activate the target site by setting s{tu + [AtJ) := 1, 
where [.J denotes truncation to an integer, and go back to step (ii). 

(v) Increment tu and proceed with the next time step. 

This particular variant of the TSP is controlled by the spreading probability p G [0, 1] 
and the Levy exponent k, G (0, 1). As can be verified easily, the assignment At := y~^/'^ 
generates a probability distribution with a lower cutoff at At = 1 which reproduces 
the asymptotic decay postulated in Eq. ([1]). Moreover, the model is defined in such a 
way that the number of attempted activations n(t„) discussed in the previous section 
occurs with probability — p), hence n = p/{l—p). Note that in this update scheme 
repeated activations of the same target site have no effect. 

The above algorithm can be easily implemented on a computer. For example, if 
rndO returns a random number drawn from a flat distribution between and 1, a 
minimal C-code for this update procedure would read as follows: 

const int Tmax=10000; // maximal cutoff time; 

const double kappa=0.5, p=0. 574262; // control parameters; 

int s [Tmax] ; // the timeline; 

for (int t=0; t<Tmax; ++t) s[t]=0; // clear timeline; 

s[0]=l; // place initial seed; 

for (int tu=0; tu<Tmax; ++tu) // execute update loop 

if (s[tu]==l) // over all active sites 

while (rnd()<p) { // repeatedly with prob. p; 

double dt = pow(rnd() , -1/kappa) ; // generate a time interval 
if (tu+dt<Tmax) s [tu+f loor(dt)] =1 ; // and activate target site. 
} 

Depending on the initial state, the algorithm can be accelerated significantly by storing 
the active sites in a dynamically generated list instead of using a static array. The 
structure of such an optimized code is outlined in the appendix at the end of this paper. 



2.2. Phenomenological properties 

Starting with a single seed of activity at t = and averaging over many independent 
runs we measured the probability g{t) to find an active site at time t. Varying the 
parameter p we observe the following phenomenological behavior (see Fig. [3]): 



6 




t 



Figure 3. Decay of the density of active sites in the TSP for k = 0.5 and various 
values of p, as listed in the legend. A qualitatively similar behavior is observed for all 
values of < K < 1. 

• For small values of p the density of active sites decays as g{t) ~ t^^^^'^K This power- 
law decay is a direct consequence of the Levy distribution ([T]) and characterizes the 
subcritical phase of the TSP. 



• For large values of p the density first decreases until it reaches a minimum, then 
increases again until it saturates at a stationary value. 

• At a well-defined critical threshold p = Pc{i^) the density decays algebraically but 
much slower as in the inactive phase. The observed power law 

Q{t) ~ (2) 

is very clean and the critical exponent a is found to vary continuously with n. It 
turns out that the numerical estimates are in excellent agreement with the analytical 
result a = 1 — ft, as will be discussed in the following section. 

For localized initial configurations such as a single seed of activity at t = the process 
may terminate when no active sites to be updated are left. In this case the last active 
site on the timeline defines the time where this particular run ends. 

Likewise one can estimate the probability Ps{t) that the process starting with a 
single seed survives at least until time t, i.e., it produces at least one active site at time 
t' > t. Averaging over many runs at criticality the survival probability seems to decay 
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algebraically as 



Ps{t) ~ t 



-5 



(3) 



although the observed scaling is less clean in this case. The survival exponent 5 is 
significantly smaller than a. For example, for k = 0.5 one finds 6 = 0.165(3). The 
unusual decay of the survival propability will be discussed in details in the last section 
of this paper. 

Continuous phase transitions into absorbing states are generically characterized by 
four independent critical exponents /3, i'^, i^y. The first two exponents are related to 
the order parameter and its conjugate field while the latter describe how the spatial and 
the temporal correlation lengths diverge as the critical point is approached. Interpreting 
the TSP as the zero-dimensional limit of directed percolation with temporal Levy fiights, 
it has no spatial degrees of freedom and hence the exponent z/^ does no longer exist. 

Moreover the question arises whether the so-called rapidity reversal symmetry, 
which in the case of DP forces the exponents (3 and (3' to be identical [19], still holds 
for the TSP. As shown in Ref. [3] this is indeed the case, despite the unusual value of 
the survival exponent 6. The purpose of this study is to find out how the exponents a, 
z/||, and 6 are related to k. 

3. Field theory 

In this section we describe a field-theoretic approach for the TSP which allows one to 
compute the critical exponents perturbatively by a loop expansion. This field theory was 
first introduced by Deloubriere and van Wijland, somewhat hidden in an appendix of 
Ref. [3] . Here we rederive their results independently, calculating the critical exponents 
to one-loop order and confirming their results. 

3.1. Langevin equation 

Before discussing the Langevin equation for the TSP, let us first to recall the well-known 
Langevin for ordinary directed percolation [20] which describes the temporal evolution of 
the coarse-grained density of active sites g{x, t) in a DP process in d spatial dimensions: 



Here the first term on the r.h.s. accounts for offspring production and spontaneous 
removal of particles. This means that the parameter a is related to the percolation 
probability and has to be tuned to a certain critical value at the transition. The second 
non-linear term prevents the density from diverging and refiects the fact that repeated 
activations of the same site have no effect. The third term describes nearest-neighbor 
diffusion in space while ^ (x, t) denotes a white Gaussian noise which accounts for the 



dtg{x,t) = ag{x,t) — bg'^{x,t) + DV^g{x,t) + ^{x,t) . 



(4) 
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fluctuations of the coarse-grained density g{x, t) caused by the stochastic nature of the 
dynamical rules. Since the intensity of such fluctuations depends on the density of active 
sites, the central limit theorem implies that the correlations of the noise are given by 

, t')) = cgix, t)5\x - x')5{t - t') . (5) 

In order to find a Langevin equation for the TSP, the DP Langevin equation given 
above has to be modified in two ways. On the one hand there are no spatial degrees of 
freedom in the present case, meaning that the spatial argument x as well as the diffusion 
term DV^g have to be omitted. On the other hand, the non-Markovian dynamics by 
means of directed Levy flights has to be incorporated. As shown in previous studies (see 
e.g. [18] and references therein), this can be done by replacing the temporal derivative 
dt on the l.h.s. by a so-called fractional derivative defined by [2] 

dt Qit) = ^ dt' t'-'-^m - git - t')] , (6) 

where k G [0, 1] is the control exponent introduced in Eq. ([T]) and A/||(/«) = — r(— k) is a 
normalization constant. The effect of this operator is to transfer activity located at time 
t — t' over a temporal distance t' > to the destination t at a rate proportional to t'~^~'^. 
Therefore, the fractional derivative generates directed Levy flights according to the 
distribution ([1]) controlled by the parameter k, just in the same way as an ordinary 
derivative dt generates a local translation in time. 

Technically the easiest way to handle the fractional derivative is to consider its 
action in Fourier space where it brings down a factor {—iuj)'^ in front of the exponential: 

d^e-"^' = {-iufe-''^' for </€<!. (7) 

In this representation the directed character of the Levy flights, which is needed to 
ensure causality of the temporal evolution, is reflected by the fact that {—iooY is not 
invariant under the replacement — cu. 

At first glance, the Fourier representation ([7]) suggests that for k = 1 the fractional 
derivative d\ should act in the same way as the ordinary derivative dt- However, it is 
important to note that for k = \ Eq. ([71) is no longer valid. In fact, d\ is a non-local 
operator while dt is not. 

With these two modifications the non-Markovian Langevin equation for the TSP 
reads 

~dtg{t) = ag{t)~hg(tf^i{t), (8) 
where ^ is a density-dependent noise with the correlations 

mm) = og{t)b{t-t!). (9) 
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3.2. Partition sum and field-theoretic action 



The partition sum Z for the present model is defined as the sum over all configurations 
and all realizations of randomness which obey the dynamical rules. In the continuum 
limit this corresponds to the functional integration over all configurations of the field 
Q{t) and all realizations of the noise ^{t) weighted according to the correlations ([9]) which 
obey the Langevin equation ([8]). Formally this may be written as 



Z oc Dq D^P[^] 6 TdtQ{tY - aQ{t) + hQ^{t) - i{t) 



(10) 

where the Langevin equation appears as the argument of a Dirac-5 functional and 

PK|o.exp(-/;jd«gl) (11) 

is the functional weight of the Gaussian noise. For later convenience we introduced an 
additional coefficient r in front of the fractional derivative which fixes the overall time 
scale of the temporal evolution. 

Following standard methods described in [21,22] one can integrate out the noise as 
follows. First the (5-functional is represented in Fourier space by introducing a response 
field Q{t): 



Z oc / DqDq / D^P[^] exp 



dt Q ( rdtg'^ — ag + bg"^ — ^ 



After a Wick rotation in the complex plane the density- dependent noise contribution 
can be separated: 

"+00 



Z (X / DgDg exp 



dt g ( rdtg'^ — ag + hg' 



X / D^P[i\ exp 



+00 



dtg^ 



This allows the noise to be integrated, resulting into 
Z (X DgDg exp 



dt[ g[Td^ — a]g + bgg — -g g 



(12) 



(13) 



sum 



With 0(t) := ^J2h/cg{t), (j)[t) := ^Jc/2h g[t) and g := \/2bc one obtains the partition 

Z oc I D^D4>e'~^^'^'^^ (14) 



with a field-theoretic action S = Sq + Sint consisting of a free part 

"+00 



So[4>,4>] 



dt (f){t) rd'^ - a (f){t) 
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(15) 



and an interaction part 



S-,r 



+00 



dt 4>{t) - 



(16) 



Apart from the fractional derivative and the missing spatial degrees of freedom, this 
action has exactly the same structure as in the case of ordinary directed percolation. 



3.3. Rapidity reversal symmetry 



As pointed out in [3] the well-known rapidity reversal symmetry of directed percolation 
still holds in the present case, i.e. the field-theoretic action derived above is invariant 
under the replacement 



0(t)^-0(-t), 0(t) 



(17) 



To see this it is convenient to represent the action in frequency space by introducing the 



Fourier transforms $(u;) := dte *'^*0(t) and likewise ^{uj), turning the action into 

'$(cj) (18) 

-UJi) (pi^UJi + UJ2) ■ 



So 



int 



- 

— $1 

'-00 27r 

g_ I d^i 
2 



-UJ 



+00 



{UJI + UJ2) 



-UJ2 



) 27r 2-71 

where we used Eq. ([7]). As can be seen both parts of the action are invariant under 
the replacement $(c<j) —^{—u) and $(c<j) — — $(— cj), hence the rapidity reversal 
symmetry still holds. This symmetry implies that the scaling dimensions of the fields 
and 4> have to be identical, i.e., /3 = reducing the number of independent critical 
exponents by one. 



3.4. Dimensional analysis 

Using the rapidity reversal symmetry the theory at tree level is expected to be invariant 
under the scale transformation 

t^Xt, (f)^X-''(f), 0^A-^0 (19) 

with a time dilatation factor A > and a field exponent x = the Fourier- 

transformed action this corresponds to the replacement 

uj^X-^uj, ^^X^-^^, <I^A^"^$. (20) 

It is easy to check that this scale transformation can be compensated by changing the 
coefficients as 

r^A^^-i+V, a-^A^^-^a, g^X^^'-^g. (21) 
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In order to establish scale invariance at tree level all coefficients have to be invariant. 
Firstly, the invariance of r implies that x = (l~'^)/2- Secondly, we have 2% — 1 = k > 
so that coefficient a is relevant, hence it has to be set to zero which is just the mean-field 
(MF) critical point. Finally, scale invariance at tree level requires the coefficient g to be 
irrelevant, i.e. 3% — 1 < or equivalently k < 1/3. Therefore the value 

i^c = l (22) 

plays the role of a lower critical threshold where mean-field behavior sets in, comparable 
to the upper critical dimension dc in ordinary directed percolation. For k, < Kc the model 
is expected to exhibit mean-field behavior with the critical exponents 

^MF ^ ^ ^ ^MF ^ ^-1 _ (23) 

The case k > Kc, where fiuctuation effects have to be taken into account, will be 
addressed in the following section. 



4. Renormalization group calculation 

4.I. Loop expansion 

In seed simulations the density g(t) measures the response of the system at time t to 
an activation of a single site at t' = and therefore can be interpreted as a two-point 
correlation function. In the field-theoretic framework this correlation function can be 
expressed as 

G{t-t') = mm), (24) 

where (. . .) denotes the statistical average according to the action. The density g{t) 
measured in seed simulations is expected to be asymptotically proportional to G(t), 
which allows one to compute the exponent a in Eq. ([2]). 

In order to compute the two-point function it is convenient to add external currents 
J{t) and J{t) to the partition function f[T^ : 

D<PD^ exp -5o[0, 0] - Si„t[0, 0] + / dt{J<P + J0) . (25) 

J — oo 

This allows the correlation functions to be expressed as functional derivatives 

(26) 



J=J=0 

As usual one separates the partition function into a free and an interacting part 

Z[J,J]cx exp(-5i„t[^,^]) (27) 

12 



X / D(pD(j) exp 



dt( J0 + J0) 



Integratinng the remaining Gaussian problem one is led to 



Z[J, J] oc exp( —S\ 



- 6 


6 ■ 








6J. 


j exp 
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J{-iu)Goiuj)J{uj) 



(28) 



where J{uj), J{uj) are the Fourier-transforms of the currents J{t), J{t) and 
Go{uj) ^ 



(29) 



denotes the free propagator. 

At this point it is important to note that the structure of the field theory is exactly 
the same as in the case of ordinary directed percolation (for a recent review see e.g. 
Ref. [23]). In particular, the loop expansion and the Feynman graphs are exactly the 
same. What changes is only the form of the free propagator and the absence of spatial 
degrees of freedom. Therefore, we can formally use the same loop integrals as in DP, 
simply omitting the integration over momenta and using the modified free propagator. 
For example, the one-loop expansion for the two-point correlation function is given by 



+ 00 



27^<2+" 



Gr 



0{g') (30) 



4-2. Wilsons renormalization group scheme 

For K, > 1/3, where the critical behavior of the TSP is influenced by fluctuations, 
the integrals in the loop expansion diverge in the limit — oo. In this case the bare 
continuum description is no longer meaningful, instead the discrete nature of the process 
has to be restored through the back door by regularizing the integrals. 

In the following we adopt Wilsons renormalization group (RG) scheme which turns 
out to be particularly suitable for the present problem. In this approach the UV 
divergences are regularized by introducing a cutoff in momentum space. In the TSP, 
where momenta are absent, we introduce a cutoff Q in the frequencies instead, i.e., the 
integration range of the loop integrals is restricted to G [—^, +^]- 

Let us now consider an infinitesimal scale transformation with A = 1 — e. Without a 
cutoff the coefficients would change according to Eq. ( 12T1) . However, the cutoff frequency 
has to be rescaled as well by — (1 + e)Q, modifying the value of the integrals in the 
loop expansion. To take this additional change into account the integrals are evaluated 
within the frequency shell |a;| G [Q, (1 -|- a process called shell integration. Finally 
the resulting contributions are expanded in to lowest order in absorbed in the coefficients 
by adding suitable terms Lr, La, and Lg on the r.h.s. of the RG equation: 

delnr = 1 — K — 2x — Lr 
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djna = l- 2x- La (31) 
djng = l- 3x- Lg 

The loop corrections L^, La, and Lg depend on the parameters r, a, g as well a the cutoff 
Q and will be computed below. 

4.3. Renormalization hypothesis of the Levy operator 

In preceding studies of models with long-range interactions by Levy flights it turned out 
that fleld-theoretic loop corrections do not renormalize the fractional derivative itself, 
instead they always renormalize the corresponding short-range operator. Technically 
this can be traced back to the fact that loop expansions always yield power series in k 
and 00 with integral powers which correspond to ordinary derivatives in real space. This 
observation implies that fractional derivatives are modifled under scale transformations 
exclusively by their scaling part, giving rise to an exact scaling relation among critical 
exponents. Assuming the same to be true in the present case, this means that L^ = 
even in the fluctuation-dominated regime k > Kc, implying the scaling relation 

X = f = ^. (32) 

II 

This scaling relation is expected to hold exactly over the full range < k < 1 to all 
orders of perturbation theory. As a direct consequence, the density in seed simulations 
Q{t), which is proportional to G(t) = {(f){t)4>{0)) , is predicted to decay as 

g{t) ~ t-^'^/^ii ~ (33) 
for any < k < 1, i.e. the decay exponent in Eq. ([2]) is given by 

a = l-K. (34) 
In numerical simulations (see below) this prediction is conflrmed with high precision. 

4.4- Analysis of the renormalization group flow 

In the second step of Wilsons RG scheme, the so-called shell integration, the propagator 
and the vertex coefficient change to one-loop order by 

G^\u) G^\u) + ^ 1^ ^Go(| + ^0^0(1 - uj') (35) 



and 



g ^ g - 2g^ I^^Glico)Goi-co) . (36) 
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Here Gq^u) = {j^—iujY — is the free propagator and ' >' denotes integration over 
the frequency shell \uj\ G (1 + e)VL\. Integrating and expanding to lowest order one 
obtains 



—Gl{uj')G^{-J) = ^2 ^ ^2^2. _ 2aT^- cos(f )) 



7r(a2 + r2^]2« - 2arfi'^ cos(^)) 
-e^{a-rQr cos(^)) 



Therefore, the loop corrections are given by 

9^^ 



" ~ 27ia (a2 + r^Q^- - 2arf]- cos (f )) 
^ 2^2j^(a-rr]'^cos(^)) 

" 7r(a2 + r2fi2«-2arr]''cos(f ))^ 

Together with the scaling relation ([32]) the RG equations read 

9,lnr = 

d^ln a = K, — La 
3k - 1 



+ 0{u^) 
(37) 

(38) 
(39) 



(40) 



djng 



2 

Their non-trivial fixed point is given by 

3fi: - 1 



l: = k. 



' 2 

or, in terms of the original paramters, by 



(41) 



kVL" (2{7k - 1) cos (^) + \/2^cos(7rK)(l - 7k)^ + (14 - 17k)k - 1 



(9 



*\2 



2^3o3k-1 



(llK- 1)3 



COS 



22k- 2 



(42) 



+ \/2cos(7rK)v/cos(7rK)(l - 7k)2 + (14 - 17k)k - 1(1 - 7k) 
+ (k(73k - 22) + 1) cos 

+ 4V2KVcos(7rK)(l - 7k)2 + (14 - 17k)k - 1 

In order to determine the exponent i/y let us consider the RG flow in the vicinity of 
this fixed point. The linearized flow field is given by the matrix 



•3k-1 



M 



daa{K - La) daa{ 



Lr. 



(43) 



a=a*,g=g* 



and can be computed explicitly 
M = 



{1-k)/4 f{K) 
—2k, 1 — 3k 



(44) 
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where 

/(«:) = 2kcsc2(-)- — -- + 2 (45) 
- ^-^ Vcos(7r/t)(l - 7/t)2 + (14 - YJk)k - 1 cot j esc j . 
The eigenvalues of this 2x2 matrix can be computed exphcitly. Defining the distance 

e = K — Kf. = K (46) 

3 

which plays a similar role as the dimensional difference e = — d \\y field theories 
with spatial degrees of freedom, and expanding the eigenvalues to lowest order in e one 
obtains 

Ai = -35 + 0{^) , A2 = i + I + 0{e') . (47) 

The first eigenvalues is always negative while the positive one describes how the critical 
parameter vanishes under rescaling. Therefore, one can identify the second eigenvalue 
with i^ii = A2 ^, leading to the main result 

i^\\= 3-^-8 + 0{e') . (48) 

Together with the exact scaling relation (!32l) this implies 

P = l-le + 0{e') . (49) 

These findings are in full agreement with the results by Deloubriere and van Wijland in 
the Appendix of Ref . [3] . Their calculation involves a dimension-like parameter d which 
is related to the control parameter k in our work by 0? = 2 — 2k. 



5. Numerical results 



We performed extensive numerical simulations with a code based on dynamical lists as 
described in the appendix. First we determined the critical parameter pc and estimated 
the exponent a for various values of k. It turned out that the scaling relation a = 1 — n 
(see Eq. ( l34l) ) is obeyed with at least three digits accuracy in the range 0.3 < k < 0.8 
so that we have no doubt that this scaling relation is correct. That is why we decided 
to consider this scaling relation as given and to use it for a precise determination of 
the critical point. The results are listed in Table [H As can be seen, the estimates are 
less accurate for very small values of k, where the Levy fiights become extremly long- 
ranged, as well as in the limit k — > 1, where the particle densities are so high so that 
the list-based algorithm is no longer efficient. 

The same type of simulations was used to measure the survival probability -Ps(t) 
at the critical point. As shown in the right panel of Fig. HI the power laws are less 
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Figure 4. Decay of the density of active sites g{t) (left) and the survival probability 
Ps{t) (right) in seed simulations of the TSP at criticality for various values of k averaged 
over 4 x 10^ up to 1.4 x 10^ runs. The density of active sites shows a very clean power 
law g{t) ^ t~'--^~'^\ The survival exponent 6 is estimated by averaging over the last 
four decades in the right panel. 
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Table 1. Numerical estimates of the percolation threshold Pc and the critical 
exponents , 6 for various values of k. Entries without error bars are based on 
analytical arguments. The estimates for S have to be compared with the predicted 
values Sc according to the conjecture in Eq. (j53p truncated to three digits. The values 
for P in the last column were computed from the previous exponents using the scaling 
relation (3 = 6i>\\ . 



clean in this case. A conservative analysis leads to the estimates reported in Table [TJ 
Apart from a slight systematic deviation inside the error bars, these estimates are in 
good agreement with the predicted value Sc which will be derived below in Eq. (!53|) . 

In order to determine the exponent uii we performed extensive off-critical 
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Figure 5. Off-critical simulations, here for k = 0.4 with p — ranging from ±0.00016 
to ±0.01024. The measured curves (left) can be collapsed convincingly (right), allowing 
one to estimate the exponents a and i^n . 
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Figure 6. Numerically measured critical exponents (black) for various values of k 
according to Table[T]compared with the mean- field prediction (blue), the field-theoretic 
one-loop approximation (red) and the conjectured formula for 6 (green), see Eq. (|53p . 



simulations. A typical data set is shown in the left panel of Fig. [51 Assuming the usual 
scaling form g{t,p — pc) — t~°'R{t{p — Pc)^") and using the scaling relation a = 1 — k 
the curves should collapse if t^g is plotted against t{p — Pc)"^^. As demonstrated in the 
right panel of Fig. [5], one obtains a very clean collapse below and above criticality. This 
allowed us to estimate the exponent z/y , as listed in Table [H 

As a visual summary the estimates for the critical exponents are plotted as functions 
of K, in Fig. [6] and compared with the analytical predictions. As expected, the field- 
theoretic one-loop expansion is tangent to the data in k, = Kc = 1/3, as indicated by 
the red lines. As can be seen, our numerical simulations fully support the field-theoretic 
results. They would even allow us to verify the results of a future two-loop calculation 
by fitting a parabola. Obviously the two-loop corrections should have the opposite sign. 
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6. Related critical properties 



In this section we discuss various issues such as the conjugate field and the role of 
different initial conditions. Based on these arguments we arrive at a conjecture that 
allows us to express the survival exponent 6 in terms of an exact scaling relation. 

6.1. Conjugate field 

In equilibrium critical phenomena an order parameter is always associated with a 
conjugate field h that, when applied externally, causes a response of the order parameter. 
For example, in the Ising model h is just an external magnetic field. The same 
applies to non-equilibrium phenomena as DP, where the external field corresponds to a 
spontaneous creation of activity at rate h. 

In the mean field regime of DP a constant external field causes an asymptotically 
stationary response pstat ~ Vh- In the non-trivial regime, where fiuctuation effects are 
relevant, the situation is different. Here a constant field h is known to cause a response 
which scales as where a = i/y + du± — (3' . In the present case, where d = Q and 

(3' = j3, this would imply that a = — (3. Together with the scaling relation (|32l) one 
arrives at 



This power law is in good agreement with numerical results (not shown here). 
6.2. Fully occupied initial state for t < 

So far we considered simulations starting with a single active seed at t = 0. However, as 
already discussed in the Introduction, the initial condition is not determined by the state 
of the site at t = alone, instead the non-Markovian dynamics requires us to specify 
the configuration of all sites along the entire time line, in principle even including all 
sites at negative times t < 0. As an example let us consider the intial configuration 



where all sites at negative times are initially occupied. This initial density is shown as a 
black step function in the left panel of Fig. [71 For this initial configuration the process 
has to start at t^ = — oo. Upon reaching t„ = the update rule will have activated 
many sites at t > 0. These active sites are randomly distributed and uncorrelated with 
an expectation value decaying asymptotically as 




for K < i 
/iT+f for K > i . 



(50) 




(51) 



g{t) = {s{t)) ~ t 



— K 



{tu = 0) 



(52) 
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Figure 7. Initial configuration where all sites at negative times are active. The left 
panel shows the initial configuration (black), the density profile in the moment when the 
update rule reaches tu = (blue) and the final density profile in the limit ^ oo (red). 
The right panel shows the same data in a double logarithmic representation. 

as indicated by the blue curve in Fig. [3 

Subsequently, as the update rule advances from t = to t — > cxo, even more sites 
will be activated along the timeline. The resulting density of active sites is shown as 
red curve in the figure. As can be seen, the density seems to decay asymptotically 
as g{t) ~ t~^/^ for k = 0.5. Obviously, this density is much larger than the density 
produced by a single seed (where one would obtain g{t) ~ t~^/^). 

Moreover, the observed exponent 1/6 coincides with the survival exponent in seed 
simulations. This relationship is valid for all k and is a direct consequence of the time 
reversal symmetry discussed in Sect. 13.31 On a quahtative level it can be explained as 
follows. As sketched in Fig. [8] a particular run in seed simulations is said to survive (at 
least) up to time t if there is a connected path of subsequent Levy flights from the seed 
to some site at t' > t. The lower part of the figure shows the same but horizontally 
reflected realization of Levy flights. In such a time-reversed situation survival translates 
into the condition that a site at time t is connected backwards in time to some site at 
t' < 0, or equivalently, that a site at time t is activated by an initial configuration where 
all sites at negative times are active. 

We therefore conclude that the survival probability in seed simulations decays 
asymptotically in the same way as the particle density in simulations starting with 
fully occupied lattice at negative times. This unusual relationship will be used in the 
following to determine the value of 5. 

6. 3. The survival exponent 5 

We now propose a conjecture for an exact scaling relation that determines the value 
of the survival exponent. The conjecture is based on the findings of the previous 
subsections and works as follows. Let us consider again the TSP starting with the initial 
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configuration flSip where all sites at negative times are active. As discussed before and 
demonstrated in Fig. [7] the process evolves in two steps: 



(i) During the updates from tu = —oo to = 0, the process produces an uncorrelated 
random distribution of active sites at t > which on average decays as t~'^. 

(ii) This random distribution can be interpreted as a time-dependent external field 
h{t) ~ t~'^, which causes the process to create additional active sites as a response 
during the updates from = 1 to = cxd. 

We now assume that this external field h{t) decays so slowly that the response of the 
process changes adiabatically as if h was stationary. Numerical experiments support the 
validity of this assumption over the full range of k. This would lead to the conjecture that 
the final density of active sites (the blue curve in Fig. [7]) decays as p{t) ~ h^^'^ ~ t~'^^^"', 
hence the survival exponent should be given hj S = K(3/a. Inserting Eq. (150|) this would 
mean that 

f ^ for K < i 
H^) = { \ (53) 
I TT^ for > i 

As can be seen in Tableland in the central panel of Fig. [3], this prediction is in excellent 
agreement with the numerical estimates for the survival exponent 5. This observation 
raises the hope that the formula (|53|) may qualify as an exact scaling relation. 




Figure 8. Illustration of the time reversal symmetry. The upper part (a) shows a 
particular realization of Levy flights in a simulation starting with a single seed. At time 
t' the process is still surviving because there is at least one connected path (bold lines) 
that activates a site at t > t'. The lower part (b) shows the time-reversed configuration 
of Levy flights which has the same statistical weight. Here the site at time t' becomes 
active if there is at least one connected path to a site at t < 0. Therefore, a fully 
activated initial state at t < (indicated as a red bar) will lead to a density profile 
that decays in the same way as the survival probability. 
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7. Summary 



In this paper we have investigated a spreading process on a timehne which exhibits a 
continuous phase transition out of equihbrium. Studying a particular model of the TSP, 
where active sites activate other sites in the future by means of temporal Levy flights, 
we estimated the critical exponenets as functions of the Levy exponent k (see Tabled]). 
Moreover, we confirmed earlier field-theoretical results obtained by Deloubriere and van 
Wijland [3] by an independent renormalization group calculation. The combination of 
field-theoretical and numerical methods lead to the following main results: 

(i) For < K < 1/3 the model exhibits mean-field behavior while for 1/3 < k < 1 
fluctuation effects (loop corrections) have to be taken into account. Here the lower 
critical value Kc = 1/3 plays a similar role as the upper critical dimension dc in 
other universality classes of continuous phase transitions. 

The TSP can be interpreted as a zero- dimensional limit of directed percolation with 
temporal Levy flights. 

im) The phenomenological scaling theory for absorbing phase transitions is still valid 
for the TSP. As there are no spatial degrees of freedom, the exponent z/^ no longer 
exists, meaning that the transition is characterized by three exponents i^y. 

liv) The time reversal symmetry, which is the essential symmetry of DP, still holds for 
the TSP and implies that f3 = [3' . 

As usual in such problems, the critical exponents vary continuously with k. 

I VI) Since the Levy operator does not renormalize itself, one obtains the exact scaling 
relation (3/v\\ = (1 — k)/2. As a consquence, in seed simualtions the density of 
active sites decay as Q{t) ~ with a = 2/3/z/|| = 1 — k. 

ivn) In the fluctuation-dominated regime a field-theoretic renormalization group calcu- 
lation to one- loop order leads to the approximation /? ^ 1 — |(k — 1/3), see Eq. (H9|) . 

(vm) Because of the non-Markovian dynamics the initial configuration requires to specify 
the state of all sites, even of those at t < 0. 

I IX) The time reversal symmetry implies that the survival probability Ps(t) ~ t^^ decays 
in the same way as the density of active sites generated by a process starting with 
a configuration where all sites at negative times are active. 
[X) Based on this observation we have conjectured that the survival exponent 6 is given 
by Eq. ( !53l) . This conjecture is in good agreement with the numerical results, see 
Fig. El 

Due to the absence of spatial degrees of freedom the fiel theory for the TSP is par- 
ticularly simple. It would be interesting to perform a two-loop RG calculation and to 
compare the results with the present numerical data. 
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Appendix: Optimized code using dynamical container classes 

In this appendix we demonstrate how the TSP can be simulated efficiently. For 
simplicity we focus on the update sequence using dynamical lists. 

Most advanced programming languages provide libraries for standardized 
dynamical container classes such as sets, lists and maps. An example is the standard 
template library (STL) which became part of most C++ environements (see e.g. [24]). 
Since the configuration of the TSP can be characterized by an ordered set of integer 
numbers marking all active sites, the suitable template class to store a configuration is 
a so-called 'set' of integers. This container class provides various functions, of which 
only four are important in the present case: 

• empty: returns true if the set is empty. 

• insert: adds a new element provided that it does not yet exist. 

• begin: returns a pointer to the first element in the set. 

• erase: removes an element from the set. 

The optimized algorithm works as follows. First one includes the relevant part of the 
STL and creates an instance of a set of integers: 

#include <set> 
set<long int> S; 

Upon creation this set is initially empty. Note that we have used the type long int as 
template argument in order to avoid possible integer overfiows caused by extremly long 
Levy fiights. Next, one has to specify the initial state. For a single seed at t = this 
can be done by 

S. insert (0) ; 

Then the update loop is executed as long as the set contains elements: 
while (not S.emptyO) tu=update (p, kappa, tmax) ; 

In the body of this loop the function 



23 



long int update (double p, double kappa, long int tmax) 
{ 

long int tu = * (S .beginO ) ; // get first element 

while (rnd()<p) // With prob. p repeat: 

{ // compute target site 

long int tnew = tu + (long int) pow(rnd() , -1/kappa) ; 

if (tnew<tmax) S . insert (tnew) ; // activate target site; 

} 

S.erase(S.begin()) ; // remove first element; 

return tu; // return update time; 

} 

performs an update and returns the time tu of the updated site. In the example shown 
above the updated site is removed from the set in order to minimize memory consump- 
tion. Note that in a set each element can appear only once; the attempt to insert an 
already existing element has no effect. In the present case this ensures that repeated 
activation of the same lattice site has no effect. 
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